We have developed general matrix-based techniques for fitting large spectroscopic data sets to complex kinetic and other models. These techniques were applied to the detailed analysis of sets of time-resolved spectroscopic data for hemoglobin measured over a range of fractional ligand dissociation, using kinetic formulations of the two-state allosteric model extended to incorporate geminate ligand rebinding and tertiary conformational changes. The first variant of the model requires three spectroscopically distinguishable tertiary states of unliganded hemoglobin subunits, two in the R quaternary structure and one in the T quaternary structure. We have refined the kinetic treatment to allow for a non-exponential form for the tertiary conformational relaxation, to which the rate of geminate rebinding is coupled; linear free energy relations connect the quaternary transition rates for tetramers with different numbers of bound ligands and/or different configurations of subunit tertiary states. This model describes the measured ligand rebinding and conformational kinetics very well and predicts a difference spectrum between the unliganded T and R quaternary states very similar to that which has been obtained elsewhere. In addition, the resulting kinetic parameters predict equilibrium properties which are consistent with the known properties. We are also studying a "tertiary two-state" model in which individual subunits in both the R and T quaternary structures exist in the same two tertiary conformations, r and t, which are identified with the high and low ligand-affinity states, respectively, of the molecule. The r conformation is stabilized by ligand binding and the R quaternary structure, and the t conformation is stabilized by ligand dissociation and the T quaternary structure. We are using specialized techniques for solving very large stiff systems of differential equations to implement a complete mathematical description of this model incorporating more than 1200 kinetic species. We are also using such "sparse-matrix" techniques to solve extremely large kinetic models (2/15 = 32768 states) describing the folding of a 16-residue $ hairpin-forming peptide in solution.